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Chapter  1 


Introduction 

1.1  Background  and  Motivation 

International  nuclear  testing  threshold  treaties,  under  effect  in  some  form  or  another  since  1963, 
have  rallied  scientists  to  devise  ways  of  monitoring  nuclear  weapons  tests.  The  UN  General  Assembly’s 
adoption  of  a  Comprehensive  Tbst  Ban  Treaty  in  1996  heightens  the  need  to  distinguish  nuclear  events 
not  only  from  earthquakes,  but  also  from  man-made  mining  explosions.  The  treaty  requires  a  monitoring 
threshold  potentially  as  low  as  magnitude  ms=2.5  (Stump  et  al.,  1994),  which  translates  to  a  yield  of 
roughly  16  T  of  explosive  (OTA,  1988).  The  CTBT  is  even  more  stringent  than  the  former  international 
agreement,  the  Threshold  Test  Ban  Treaty,  which  already  imposed  a  celling  of  150kT  (OTA,  1988;  personal 
communication  from  Dr.  Michael  Hedlinand  Professor  John  Orcutt).  In  theU.S„  the  number  of  man-made 
mining  explosions  with  yield  greater  than  50  T  is  about  10,000  per  day,  about  one  of  which  is  greater  than 
200  T  (Stump  et  aJ„  1994),  giving  ample  justification  for  a  discriminant  that  sifts  earthquakes  and  nuclear 
explosions  from  mining  explosions  that  may  occur  on  the  same  magnitude  scale.  What  is  more,  there 
are  various  ways  to  evade  the  detection  of  a  nuclear  test,  such  as  detonating  the  bomb  in  an  underground 
cavity  or  during  an  earthquake  or  certain  types  of  mining  explosions  (OTA,  1988;  personal  communication 
from  Dr.  Michael  Hedlin).  That  evasion  tactics  can  effectively  reduce  the  perceived  magnitude  by  half 
(OTA,  1988)  underscores  the  need  for  high-gain,  high-confidence  identification  methods.  If  testing  is  to  be 
monitored  globally  and  systematically,  it  is  desirable  to  differentiate  events  independently  of  local  testing  or 
mining  practices  and  local  geological  characteristics.  Scientists  therefore  seek  a  seismic  discriminant  both 
highly  accurate  and  highly  reproducible,  regardless  of  specific  source/receiver  setting. 

The  most  prevalent  form  of  mining  explosion  is  what  is  known  as  a  ripple-fired  explosion  (Lange- 
fors  and  Kihlstrom,  1978),  A  ripple-fired  explosion  consists  of  a  grid  of  explosive  shots  that  are  detonated 
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sequentially.  Michael  Hedlin  and  colleagues  have  noted  that,  almost  across  the  board,  frequency  tuod* 
illations  of  high  spectral  amplitude  and  long  time-duration  result  from  ripple-firing  and  distinguish  these 
mining  events  from  single-shot  explosions,  These  scientists  have  also  devised  a  theoretical  model  and  a 
time  and  frequency-dependent  multi  taper  transform  that  successfully  reproduce  and  display  the  modulation 
(Hedlin  et  al„  1989;  1990).  This  discriminant  has  the  advantage  that  it  appears  independently  of  specific 
source/receiver  setting.  However,  there  are  a  few  cases  where  it  may  not  appear  (Hedlin,  1997). 

I  have  sought  to  further  the  analysis  by  replacing  muHitapers  with  the  more-refined  wavelets, 
by  using  phase  differences  between  data  components  to  observe  signal  polarization,  and  by  extending  the 
analysis  to  lower  frequencies  (to  between  .05  and  1  hz).  Until  recently,  researchers  In  the  seismic  verification 
field  have,  at  be$t.filtered  time  series  with  muHitapers,  and  for  many  applications  this  is  adequate.  Taper 
transforms  provide  robust  spectral  estimates  but  a  limited  resolution  prescription.  Like  multitapers,  wavelets 
con  be  mode  orthogonal.  But  unlike  muHitapers,  wavelets  provide  a  flexible  trade-off  between  time  and 
frequency  resolution  across  the  time  and  frequency  plane,  a  feature  well- suited  to  the  wide  range  of  energetic 
frequency  content  and  the  closely-spaced  arrivals  present  in  seismic  signals.  By  guaranteeing  the  optimal 
tradeoff  between  time  and  frequency  localization,  wavelets  can  resolve  simultaneously  in  time  and  frequency 
two  pulses  both  very  close  in  time  and  very  dose  in  frequency,  whereas  the  taper  will  be  limited  to  resolving 
the  pulses  either  in  time  or  in  frequency  (Daubechies,  1987).  Therefore,  a  wavelet  can  be  adjusted  to  achieve 
a  superior  time  and  frequency  resolution  tradeoff  than  a  mullitaper  for  a  chosen  frequency  band  of  interest 
(the  tradeoff  will  be  worse  than  a  taper's  elsewhere,  but  better  for  the  chosen  bandwidth).  This  characteristic 
is  potentially  useful  fur  discriminating  spectral  features  In  the  time  and  frequency  plane,  particularly  when 
dealing  with  regional  distances  and  surface  wave  spectra,  often  beset  with  closely  spaced  pulses  and  a  great 
deal  of  scattered  energy.  Furthermore,  since  complex  wavelets  retain  the  time  and  frequency  evolution  of 
the  relative  phase  between  data  recorded  along  different  directions  in  the  earth,  we  can  reconstruct  the  time 
and  frequency  evolution  of  the  signal's  polarization.  We  have  shilled  the  analysis  to  a  lower  frequency 
range  with  the  expectation  that  some  form  of  time-independent  modulation  will  appear  even  in  those  cases 
where  it  fails  to  appear  at  higher  frequencies  and  even  when  higher-frequency  components  have  attenuated. 
Banding  in  this  range  stems  from  finite  (as  opposed  to  infinitesimal)  source  duration.  The  pattern  of  signal 
polarization  across  the  time  and  frequency  plane  provides  a  further  constraint  upon  phase  identifications  and 
upon  identification  of  source  mechanisms  underlying  spectral  features  observed  in  the  time  and  frequency 
plane.  This  may  be  particularly  important  for  regional  and  low-  to  mid-frequency  data,  again  fbr  the  reasons 
mentioned  above. 
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1.2  Thesis  organization 

The  main  objective  of  this  thesis  is  to  evaluate  results  regarding  potential  distinctions  between 
nuclear  and  mining  explosion  spectra,  and  therefore  shall  not  attempt  to  cover  distinction*  between  nuclear 
explosion  and  earthquake  spectra.  Since  seismic  measurements  of  single-explosion  calibration  shots  at 
mines  ore  easier  to  obtain  than  nuclear  explosions,  and  since  ripple- tired  mining  blasts  are  the  most  common 
form  of  mining  operation,  this  problem  is  here  recast  as  a  comparison  between  single  calibration  shots 
and  ripple-fired  quarry  shots.  1  shall  therefore  begin  by  describing  ripple-flring-what  it  is,  and  how  it  is 
implemented.  It  is  then  necessary  to  cover  the  dual  time-frequency  decomposition  known  as  a  sonogram, 
since  this  is  the  rudiment  of  our  tool,  and,  in  so  doing,  I  shall  discuss  some  of  the  features  of  ripple-fired 
mining  blasts  that  sonograms  have  made  observable.  Up  until  now,  sonograms  have  bocn  created  with 
multitaper  filters,  and  so  I  shall  discuss  briefly  their  derivation  and  important  features.  More  importantly, 
I  will  then  attempt  to  provide  a  summary  of  (he  salient  properties  of  wavelets  and  how,  following  Lilly 
and  Park  (1995),  wc  made  use  of  a  wavclct-like  filter  derived  from  multitapers.  A  description  of  the  data 
sets  follows.  Based  on  multitaper  and  wavelet  analyses  of  a  subset  of  the  data,  a  comparison  is  then  made 
between  the  resolution  in  both  cases. 

Chapter  3  presents  the  derivation  of  the  wavelet  polarization  by  way  of  a  singular* value  decom¬ 
position  and  a  summary  of  the  results  simulated  from  synthetic  data.  It  then  introduces  samples  of  results 
on  actual  data  in  order  to  help  assess  the  tool's  validity  and  usefulness  In  Identifying  wuvetypes. 

Chapter  4  launches  a  discussion  of  the  low-frequency  analysis  and  some  of  the  most  interesting 
features  found  there,  and  some  possible  ramifications  for  seismic  verification.  Particular  attention  is  paid 
to  theoretical  underpinnings  and  evidence  for  possible  differences  in  surface-wave  dispersion  in  the  two 
explosion  types.  I  discuss  the  Iow*froqucncy  results  furNRDC  and  for  the  Wyoming  experiment  in  relation 
to  this.  In  the  end  I  try  to  summarize  what  has  been  accomplished  and  what  learned  from  this  exercise,  and 
what  remains  to  be  done.  I  also  make  suggestions  for  fruitful  future  research  along  these  lines. 
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Chapter  2 


Spectra  of  single-shot  and  ripple-fired  explosions 
2.1  Ripple-firing 

Ripple- firing  is  mining  jargon  for  the  common  mining  technique  of  fracturing  rock  with  a  cascade 
of  explosions.  Single  explosives  are  spaced  in  a  generally  shallow  array  of  rows  and  columns  on  the 
order  of  tens  of  meters  apart  and  detonated  sequentially,  ideally  at  regularly-spaced  intervals  (Figure  2.1). 
Ripple-firing  is  an  efficient  means  of  loosening  large  sections  of  rock  while  minimizing  ground  motion  and 
thrown  rock  (Dowding,  1985).  It  is  the  most  common  form  of  mining  operation  worldwide  (Langefors  and 
Kihl  strOm,  1978;  Stump  et  al,  1989). 

We  can  model  the  idealized  ripple-fired  disturbance  *(t)  by  a  finite  duration  series  of  impulses 
each  separated  by  a  period  T  and  altogether  lasting  a  time  D: 

=  (»■») 

where  u>(t),  the  basic  waveform,  is  convolved  with  the  finite  shah  function  to  simulate  repetition,  and  the 
shah  function  III  is  made  finite  by  the  boxcar  function  B,  In  real  life,  however,  the  time  interval  between 
shots  strays  from  the  extrapolated  regular  values.  For  example,  projected  and  actual  shot  times  can,  in  parts 
of  the  U.S.,  deviate  by  up  to  34%  (Hedlin  et  al.,  1989).  For  this  more  general  case,  we  may  write 

a(t)  =  ][>(t-^)  (2.2) 

i 

where  STj  gives  the  time  of  the  jth  shot  relative  to  the  time  of  the  first  shot  arid  which  can  include  any 
deviations  from  a  standard  time  offset  T. 

If  the  shots  occur  in  a  spatial  array,  as  is  often  the  case,  as  in  Figure  2.1,  we  can  incorporate  (he 
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Figure  2.1  A  cross* sectional  diagram  of  a  typical  shot  array  for  a  ripple* fired  mining 
explosion.  Credit:  Dr.  Michael  Hedlin,  IGPP*SI0. 

contributions  made  to  the  apparent  shot  interval  by  the  differences  in  distance  between  each  shot  and  the 
receiver.  Assuming  an  absence  of  shot-time  irregularities,  we  can  separate  apparent  offset  times  from  actual 
offset  times  by  using  the  expression 


ST*  —  p(Szjiin6  +  SyjcotB)  (2.3) 

where  6  is  the  source  to  receiver  azimuth,  p  is  the  horizontal  slowness,  and  Szj  and  Syj  are  the  actual  z  and 
y  spatial  offsets  of  the  shots  in  the  array  (see  Figure  2.2)  (Hedlin  et  al.,1990;  Smith,  1989). 

2.2  Spectra  of  ripple-fired  explosions 

The  Fourier  power  spectrum  reveals  time*independent  frequency  modulation.  Reamer,  et  al. 
(1992)  has  shown  experimentally  that  linear  superpositioning  of  shots  is  more  than  sufficient  to  explain 
amplitude  modulations  in  the  observations,  even  though  there  has  been  evidence  that  it  may  be  necessary 
to  include  non-linear  interactions  between  sources  in  the  elastic  (very  near*source)  region  (Reamer  et  al., 
1992).  The  long  life  of  the  modulations  stems  from  the  protraction  of  the  event:  the  panoply  of  wavetrains 
and  tbeir  scattered  waves  that  arrive  from  one  shot  continue  to  also  arrive  from  all  the  shots.  It  is  possible 
to  infer  the  intershot  times  from  the  bandwidth  of  the  lobe  spacing,  though  irregularities  in  time  offset  often 
obscure  this  relationship.  The  Fourier  transform  of  expression  2.2  gives 

X(f)  =  FT(£w(t-STj)}  (2.4) 
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and  the  power  spectrum  is 

P(f)  =  |Hr{/)|*((2  iiniirfSTj)1  +  (£cot2*fSTj)') 
'  j  i 


(2.6) 


(2.®) 


(Hedlin  et  al.,  1989;  1990).  The  factor  of  sines  and  cosines  gives  rise  to  the  observed  spectral  modulations, 
as  shown  in  Figure  2.3.  Figure  2.4  shows  the  spectrum  from  a  synthetic  of  50  shots  with  randomized  time 
offsets  of  mean  63  mi  and  variance  6.3  mi  (Hedlin  et  al.,  1989).  Figure  2.5,  on  the  other  hand,  results  from 
a  synthetic  of  39  cospatial  shots  spaced  at  regular  time  intervals  of  25  mi  (Hedlin  et  al,  1990).  Because 
in  the  latter  example  there  is  a  clean  relationship  between  offset  interval  and  lobe  spacing,  we  are  able  to 
deduce  from  the  lobe  spacing  of  40  Ax  (where  the  boxcar’s  transform,  i.e.the  sine  function,  has  its  peaks)  a 
shot  interval  of  25  ms.  or  1/40  A*-1.  From  the  widths  of  the  side  lobes  we  are  able  to  deduce  the  duration  of 
.975 1  of  the  entire  blast.  The  former  example  still  displays  a  strong  frequency  modulation  at  the  expected 
frequencies  that  are  multiples  of  1/63  m»“* .  i.e.,  16, 32  and  48  hx,  but  now  the  scatter  begins  to  obscure 
the  peaks  at  the  higher  overtones.  Changes  in  apparent  offset  time  due  to  spatial  offsets  in  an  array  only 
exacerbate  the  effects  of  the  scatter  and  can  introduce  modulation  into  the  lobes  that  makes  it  possible 
to  incorrectly  infer  the  shot  intervals.  (This  Is  demonstrated  in  the  synthetic  spectrum  In  Figure  2.6,  for 
example,  from  which  we  would  infer  a  shot  spacing  of  170  m*  «  1/6  A*"1,  but  for  which  in  fact  the  whole 
set  of  shots  lasts  384  mi «  l/(2,6)  A*"1 .  The  difference  in  travel  times  helps  reinforce  the  side  lobes.)  In 
fact  it  has  been  speculated  that  most  of  the  modulation  we  witness  in  the  spectra  are  not  in  fact  due  to  the 
ripple-firing  per  se,  but  to  the  finite  (as  opposed  to  infinitesimal)  event  duration  (Hedlin  et  al.,  1990). 


2.3  Sonograms 

In  order  to  observe  the  time-independence  of  the  modulation,  it  is  necessary  to  create  a  time  and 
frequency-dependent  spectrum,  known  as  a  sonogram.  This  is  done  by  taking  the  discrete  Fourier  transform 
of  windows  of  the  time  series  that  are  staggered  in  time,  which  serves  to  preserve  information  in  lime.  From 
the  result  we  compose  a  matrix  whose  elements  are  the  series’  windowed  power  spectrum  at  each  time 
and  frequency  point  (Hedlin  et  al.,  1989. 1990).  Since  multi-tapers  are  generated  in  the  same  manner  as  1 
will  show  below  in  the  description  of  wavelets  (Park,  1987),  I  will  omit  that  derivation  here.  An  example 
of  time-independent  banding  is  displayed  in  Figure  2.7,  a  sonogram  calculated  for  a  Soviet  quarry  blast. 
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Figure  2.2  The  geometry  for  calculating  array  spatial  offsets  as  they  affect  the  source* 
receiver  distance  function  and  hence  travel  time. 

Hedlin  et  al.  ruled  out  a  source/receiver  origin  (e.g.,  a  resonance  layer)  of  the  modulation  by  comparing 
data  from  different  source  types  measured  at  the  same  station  and  data  from  the  same  event  measured  at 
different  stations.  The  effect  was  also  observed  in  Norway  from  Norwegian  mining  blasts,  so  that  it  appears 
region-independent  (Hedlin  et  al.,  1990).  1b  help  isolate  the  fundamental  spectral  patterns  from  fluctuating 
amplitude  information,  Hedlin  et  al.  have  reduced  the  sonogram  to  a  binary  sonogram.  In  this  form,  the 
difference  between  a  highly  smoothed  (convolved  with  boxcar  of  width  2.5  bz)  and  a  less-well-smoothed 
(convolved  with  a  boxcar  of  width  1.  bz)  version  of  the  original  spectrum  is  flattened  out,,:  values  above 
and  below  the  local  average  are  assigned  values  of +1  and  -1,  respectively  (Hedlin  et  al.,  1989;  1990).  Tbe 
binary  sonogram  corresponding  to  Figure  2,7  is  shown  in  Figure  2.8. 

2.4  Filtering,  multitapers,  and  wavelets 

In  this  research  we  have  sought  to  exploit  some  of  the  advantages  wavelets  confer  upon  multi¬ 
dimensional  transforms  of  data.  A  two-dimensional  wavelet  transform  of  a  time  series  provides  improved 
tradeoff  between  time  and  frequency  resolution.  A  wavelet  transform  is  therefore  particularly  advanta¬ 
geous  in  the  analysis  of  waveforms  fraught  with  abrupt  changes  in  time  or  with  narrow,  higbly-energetic 
bandwidihs,  common  features  of  seismic  waveforms.  Unlike  the  Fourier  transform,  the  wavelet  transform 
does  not  restrict  the  tradeoff  between  time  and  frequency  localization  by  fixing  the  time  window  (Rioul 
and  VetterH,  1991;  Hlawatsch  and  Boudreaux-Bartels,  1992;  Lilly  and  Park,  1995;  ChakrabortyandOkaya, 
1995).  Instead  it  allows  inversely-varying  time  and  frequency  resolutions  as  It  traverses  the  time  series  and 
therefore  negotiates  the  best  time-frequency  resolution  for  the  lime-frequency  region  of  interest.  We  use 
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Figure  3.3  Spectrum  of  linear  superposition  ofl  to  5  subshots  at  the  same  location  In 
space  but  separated  by  100  mi  from  one  another.  This  intcrahoi  delay  yields  modulations 
snar/vl  everv  tO/iz  frc/lli:  Hr  Michael  HnHIin.  IfiPP-SFO 
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the  wavelets  championed  by  Lilly  and  fork,  which  establish  a  four-fold  advantage:  while  maintaining  the 
resolution  properties  of  classical  wavelets,  they  comprise  a  set  of  independent  estimates  that  can  be  averaged 
to  reduce  variance;  at  the  same  lime,  they  mimic  the  spectral -leakage  minimization  of  multitapers  without 
increasing  bias  (Lilly  and  Park,  1995). 

The  classical  wavelet  is  defined  to  be  one  member  of  the  set  of  self-similar,  translated,  and 
rescaled  versions  of  a  generic  mother  wavelet  whose  form  is  chosen  to  best  window  the  given  time  series 
(see  Chakraborty,  and  Okaya,  1995,  to  learn  bow  this  is  done).  If  the  mother  wavelet  in  continuous  time  is 
p(t),  then  the  discrete-lime  wavelets  used  to  convolve  the  time-series  are  the  set 

= «.«(‘-^)  m 

witht  -  mdt,m  =  ■  •  -,  0,  -  -  >,  Their  convolution  with  the  lime  series*  is  then  given  by 

*(«.&)  =  (2.8) 

m 

Here  a  is  the  scale  parameter,  b  a  time  translation  parameter  proportional  to  a,  and  e  an  amplitude  factor 
dependent  on  a.  dt  is  the  sampling  interval.  The  choice  of  a  sets  the  time  scale  of  the  decomposition  (a 
feature  that  is  mirrored  in  b'%  proportionality  to  a):  larger  scales  a  (and  hence  longer  wavelets)  canvass  larger 
time  intervals,  thus  revealing  more  global,  large-scale  features  of  the  signal,  and  vice  versa,  smaller  scales 
inspect  the  signal  more  microscopically  (Hlawatsch  and  Boudreaux-B  artels,  1992;  Rioul  and  Vetterli,  1991). 
a  is  also  proportional  to  the  center  frequency  U  of  the  analyzing  wavelet  and  inversely  proportional  to  /«, 
the  target  frequency  In  the  signal,  i.e.,  a  =  f0/fc  (Hlawatsch  and  Boudreaux-Bartels,  1992).  In  practice, 
discrete  wavelets  are  usually  generated  by  applying  a  filter  bank  to  a  time  series.  The  scale  parameter  a  is 
then  defined  as  o£, ;  a  member  of  (he  set  of  integers  (Chakraborty  and  Okaya,  1995). 

If  AT  is  the  lime  window,  then  AT  a  Mdt,  where  M  is  an  integer.  AT  therefore  satisfies  the 
uncertainty  relation  AT/„  =  (Mdt)fv  -  p,  where  the  constant  p  is  the  time-bandwidth  product.  The 
so-called  lime-bandcenter  product  p«  =  feMdt  is  also  a  constant.  Therefore,  -  p/pe  =  a  constant 
(Hlawatsch  and  Boudreaux-Bartels,  1992;  Rioul  and  Vetterli,  1991;  Lilly  and  Park,  1995).  This  shows 
that  the  width  of  the  time  window  varies  inversely  to  the  width  of  the  frequency  band,  which  in  turn 
varies  proportionally  to  /«.  Thus  we  obtain  good  frequency  resolution  at  low  to  mid  frequencies,  but  poor 
frequency  resolution  at  high  frequencies,  and  likewise,  good  time  resolution  at  high  frequencies,  poor  lime 
resolution  at  low  frequencies.  We  can  adjust  the  scale  to  accomodate  time  and  frequency  localization  of 
waves  containing  sudden,  short-lived  variation  in  time  or  containing  narrow  energetic  bandwidths  (Lilly 


TT  ‘d 


dnoao  unoao  ois 


2062  1?£S  6T9 


i0  :0T  <L66T-£T-TT 


50 


40 


frequency  (Hz) 


Pigure  2.4  Spectrum  of  synthetic  of  50  shots  of  randomized  time  offsets  of  mean  63  mt 
and  variance  6.3  m».  Credit:  Dr.  Micaie  Hedlin,  IGPP-SIO. 
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Figure  2.5  Spectrum  from  synthetic  consisting  of  39  cospatial  shots  spaced  at  regular 
time  intervals  of  25  mr,  the  Green’s  function  for  this  synthetic  is  the  single  calibration 
explosion  Cbemex  2  examined  in  this  paper.  Credit:  Hedlin,  1991. 

and  Park,  1995).  Figure  2.9  bears  witness  to  the  difference  between  the  lime-frequency  resolution  of  a 
Fourier  spectrum  and  that  of  a  wavelet  transform,  Fourier  and  wavelet  decompositions  are  both  applied  to 
the  function  S(t  - 1  j ) + f  (t  - 1  j)  +  eap(i2n7i  t)  +  egp(i2irf3t),  which  in  effect  images  the  resolving  power 
of  each  filter  in  both  the  time  and  frequency  domains.  The  top  figure  illustrates  the  Fourier  resolution,  and 
the  bottom  figure  the  wavelet  resolution.  Apparent  are  the  fixed  relationship  between  time  and  frequency 
resolution  In  the  Fourier  transform,  and  the  flexible  tradeoff  between  time  and  frequency  resolution  across 
the  time-frequency  plane  in  the  wavelet  transform.  Notice  the  wavelets*  improved  time  resolution  at  high 
frequencies  and  improved  frequency  resolution  at  low  frequencies  (Hlawatsch  and  Boudreaux-Bartels, 
1992). 

We  adopt  Lilly  and  fork’s  method  of  constructing  the  basis  wavelets  by  solving  the  eigenvalue 
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Figure  2.6  A  spectrum  demonstrating  changes  introduced  into  modulation  intervals  by 
spatial  offsets  between  shots  in  an  array.  Credit:  Hedlin,  1991. 

equation  maximizing  the  wavelets’  energy  spectrum  within  the  chosen  bandwidth.  Like  them,  we  will  also 
prescribe  (he  parameters  p,  /*,  and  p,  instead  of  the  parameters  a,  6,  and  M,  If  {wm}  are  the  desired 
wavelet  functions,  then 

M/i 

W(f)  =  £  wmtuit"  (2.9) 

ma-Sf/J+l 

where  1/1  <  jfe,  /<  is  the  target  frequency,  and  \f  -  /«)  <  We  want  to  maximize  the  ratio 

rf , , ,  _  'iMl  (210) 

This  is  equivalent  to  the  eigenvalue  equation 

(C+  -  C~)x  -  Ax  —  0  (2.11) 
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Figure  2.7  Sonogram  showing  regular,  time*independent  modulations  along  ihe  frequency 
domain;  the  event  is  the  Soviet  quarry  blast  examined  in  this  paper.  Credit:  Dr.  Michael 
Hedlin,  IGPP-SIO. 
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Figure  2.8  Binary  sonogram  of  the  same  Soviet  quany  blast  presented  in  Figure  2.7. 
Credit:  Hedlln,  1991. 
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where  the  MzM  matrix  C*  is 


±  _  sin  (2ff(/«±/w)(m-m,)<ft) 

*(m  -  m')dt  lZ1ZJ 

The  M  eigenvectors  of  this  eigenvalue  equation  are  the  M  M -component  wavelets  {wm}  (Lilly  and  Park, 
1995).  We  also  retain  signal  phase  information  by  forming  complex  wavelets  from  even  and  odd  pain  of 
wavelets  (Lilly  and  Park,  1995): 


W.f.n  +  iW'U 

V2 


(2.13) 


(see  Figure  2.10  for  some  of  the  wavelets  we  have  calculated).  The  procedure  of  constructing  wavelets  and 
convolving  with  the  data  is  repeated  for  each  time  increment  in  the  data  series  and  for  each  frequency  band. 

The  eigenvectors  comprise  a  set  of  orthogonal  wavelets  that  depend  on  the  target  frequency,  the 
bandwidth,  and  the  size  of  the  time  window.  Though  the  eigenvectors  we  have  defined  depend  on  the 
same  parameters  as  a  standard  wavelet,  they  do  not  constitute  exactly  self-similar  versions  of  a  prototype; 
however,  they  can  be  made  to  approach  self-similarity  if  p  and  pe  are  fixed  and  M  is  large.  Numerical 
spline  approximations  to  these  functions  are  in  fact  self-similar  (Lilly  and  Park,  1995).  Hence,  in  principal, 
we  generate  this  type  of  wavelet  in  the  same  way  we  create  an  optimal  multi-taper  function,  only  we 
parameterize  it  differently.  Moreover,  as  defined  above,  our  wavelets  are  not  only  time-limited  and  nearly 
band-limited,  but  also  orthogonal  and  therefore  statistically-independent  transforms,  permitting  spectral 
averaging  and  hence  reduced  variance.  We  have  not,  however,  been  forced  to  increase  the  size  of  /„  or 
of  AT  in  order  to  achieve  this  reduction  in  variance  (Lilly  and  Park,  1995).  Furthermore,  by  including 
both  negative  and  positive  frequency  intervals  in  the  energy  minimization  condition,  we  have  minimized 
the  contributions  of  overlapping—  and  beoce  non-onhogonal — sidelobes  belonging  to  the  leakier  wavelets, 
i.e.,  we  have  minimized  bias,  which  bounds  typical  multitaper  transforms  (Lilly  and  Park,  1995;  Riedel 
and  Sidorenko,  1995).  Only  a  subset  of  the  wavelets  will  guarantee  good  energy  confinement.  Since  this 
approach  includes  negative  frequency  intervals,  A  /  =  4/„ ,  and  about  the  first  L  =  A}  AT  -  ifuAT  -  4 p 
wavelet  transforms  will  have  energy  concentrated  in  the  bandwidth.  Therefore,  we  form  data  transforms 
with  the  L  wavelets  that  preserve  adequate  time-bandwidth  confinement  and  form  a  vector  from  these  L 
transforms: 


=E«eV. 


(J.14) 
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(2.16) 


We  esiimate  the  spectrum  II  by  averaging  the  inner  product  of  the  ^-component  vector  G  (Lilly  et  al): 


n  =  2 


(2.1«) 


2.5  The  seismic  data  set 

For  this  study  we  employed  two  seismic  data  sets:  one  from  a  temporary  seismic  array  set  up 
in  the  Black  Thunder  mining  region  of  eastern  Wyoming,  and  the  other  from  the  BAY  (Bayanaul)  and 
KJCL  (Karkaralinsk)  stations  of  the  NRDC  seismic  network  in  Kazakhstan  (see  Figures  2.11  and  2.12). 
The  Wyoming  data  are  three-component  recordings  of  velocity??  taken  by  surface  instruments??  at  a 
sampling  rate  of  1/100  <?#.  The  comer  frequency  is  .0083  ht.  The  local  geology  is  characterized  by  .  The 
NRDC  data  are  also  three-component  recordings  of  velocity  made  by  surface  and  borehole  instruments  at  a 
sampling  rate  of  1/250  epe,  the  surface  instruments  with  a  flat  response  down  to  1  hz  and  the  borehole  with 
a  flat  response  down  to  .2  ht.  This  area  of  Kazakhstan  is  seismically  quiet  and  tectonically  stable,  and  the 
instruments  situated  on  granitic  outcrops  to  reduce  reverberations  (Hedlin,  1991). 

The  1996  Wyoming  experiment  was  conducted  explicitly  to  gather  regional  three-component 
array  recordings  of  known  quarry  blasts  detonated  at  the  Black  Thunder  coal  mine,  This  data  set  provides 
a  generous  set  of  ripple-fired  mining  recordings  at  a  variety  of  azimuths  encircling  the  events  (see  Figure 
2.11).  Here  we  have  analyzed  the  data  taken  on  July  19  of  a  quaxiy  blast  consisting  of  a  total  of  620 
shots  arranged  in  7  rows  and  interdelayed  by  about  35  me.  The  data  from  the  BAY  and  KKL  stations 
of  the  NRDC  network  were  chosen  because  they  are  three-component  recordings  of  known  single-shot 
and  ripple-fired  mining  explosions  at  a  local  quany  and  have  a  good  signal-to-noise  ratio.  We  include  In 
our  analyses  a  known  ripple-fired  event  (day  135  of  the  1987  recordings)  and  20  ton  single-shot  mining 
calibration  explosion  known  as  Chemex  2  (day  245  of  the  1987  recordings)  so  as  to  compare  recordings  of 
single  and  ripple-fired  source  types  in  the  same  region  and  at  the  same  stations.  Though  differing  regional 
geologies  and  specific  mining  practices  present  drawbacks  for  comparisons  between  the  two  data  sets,  their 
similar  sourcc-to-receiver  range  presents  one  advantage.  We  require  three-component  data  in  all  cases  in 
order  to  perform  the  polarization  calculations,  and  we  desire  recordings  at  different  azimuths  in  order  to 
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time  step  (total  =  57  for  this  window  and  bandwidth) 


Figure  2.10  An  example  of  three  consecutive,  real  wavelets.  Note  their  alternating  even 
and  odd  symmetry  and  that  these  have  not  been  interpolated. 


observe  the  azimuthal  dependence,  if  any,  of  spectral  characteristics,  such  as  frequency  dependence  and 
time-independent  banding.  Before  any  of  (he  data  are  analyzed,  the  three  orthogonal  lime  series  are  first 
rotated  into  an  orthogonal  radial-transverse-vertical  coordinate  system  defined  by  the  ray  path  of  the  signal 
(radial  along  the  dirciion  of  the  ray  path,  transverse  perpendicular  to  the  ray  path  and  parallel  to  the  earth's 
surface,  and  vertical  perpendicular  to  the  earth's  surface).  The  analysis  at  low  frequencies  benefits  from 
first  filtering  and  decimating  tbe  data.  Both  the  NRDC  and  the  Wyoming  data  have  been  decimated  to  1 
point  in  25.  The  NRDC  data  have  been  low-pass  filtered  below  2  hi,  and  the  Wyoming  data  below  1  ht,  in 
both  cases  using  a  1  ha  lowpass  convolution  filter,  Important  geological  parameters,  such  as  local  crustal 
depths,  densities,  etc.,  are  summarized  in  Table  1.1. 


2.6  A  comparison 
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station 

H 

P 

a 

P 

km 

gm/em * 

km/  » 

km/» 

Wyoming 

39 

2.85 

6.30 

3.56 

NRDC 

41 

2.90 

6.48 

3.65 

1 . 1  Average  crustal  parameters 


Figures  2.13  and  2.14  display  the  spectral  resolution  characteristics  at  high  frequency  ranges  of 
themultitaperand  wavelet  transforms,  respectively.  The  wavelet  transforms  were  computed  using  p*  =  9.0 
and  p  =  2.5  to  favor  frequency  resolution  as  much  as  possible  over  the  wide  range  of  frequencies  under 
consideration.  It  is  important  to  note  when  observing  these  figures  that  Figure  2.14  examines  a  smaller 
subsection  of  Figure  2.13.  The  wavelet  transforms  were  computed  using  pt  -  9.0  and  p  =  2.5  to  favor 
frequency  resolution  as  much  as  possible  over  the  wide  range  of  frequencies  under  consideration  without 
compromising  excellent  lime  resolution.  The  multitaper  transforms  were  computed  using  IS  «  windows 
and  p  =  4.0.  At  these  parameter  values,  the  wavelets  drastically  improve  time  resolution  at  all  times, 
the  accuracy  being  of  course  worst  at  the  lowest  frequencies  analyzed.  The  wavelet  frequency  resolution 
may  appear  to  be  worse  than  the  multitaper  frequency  resolution,  but  I  believe  It  Is  in  fact  not  worse.  It 
only  appears  to  be  inferior  because  it  is  poor  relative  to  the  excellent  time  resolution,  whereas  the  ratio  of 
frequency  resolution  to  time  resolution  in  the  multitaper  spectrum  is  much  closer  to  one.  Therefore,  at  this 
value  of  Pc,  the  wavelets  do  not  improve  the  frequency  resolution  achieved  with  multitapera,  though  at  this 
value  of  p  they  do  improve  the  lime  resolution.  Figure  2.15  displays  the  wavelet  spectrum  of  the  same 
data  calculated  with  the  same  value  of  p  but  with  p«  =  25.0.  The  frequency  resolution  shows  considerable 
improvement  over  its  previous  values,  and  finally  appears  to  be  on  a  par  with  the  time  resolution.  However, 
such  high  values  for  pa  enforce  prohibitively  large  and  time-consuming  calculations  at  the  low  frequencies 
(below  4  h»)  that  we  examine  in  this  study.  The  choice  p*  =s  9  appears  to  be  about  the  limit  to  which  we 
can  stretch  our  demands  given  the  numerical  constraints  imposed  by  low-frequency  calculations  using  the 
Matlab  software.  It  is  interesting  that,  though  the  time-independent  banding  is  still  observed  with  improved 
time  resolution,  it  is  less  pronounced.  Multitapersmay  therefore  prove  advantageous  if  a  visual  recognition 
of  this  pattern  is  desired.  The  muititapers  give  better  frequency  resolution  at  high  frequencies,  and  wavelets 
at  low  frequencies.  Likewise,  at  the  lowest  frequencies,  the  multitaper  resolution  is  almost  on  a  par  with 
the  wavelet  resolution.  This  suggests  that  wavelets  are  advantageous  for  frequency-dependent  analyses  at 
lower  frequencies,  and  muititapers  are  perhaps  still  adequate  for  time-dependent  localization  and  separation 
at  lower  frequencies  and  that  in  this  case  a  multi-pronged  approach  may  be  optimal. 
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Chapter  3 


Wavelet  Polarization  Analysis 


3.1  Wavelet  singular- value  decomposition 


Complex-valued  wavelets  give  a  means  of  evaluating  a  wave’s  polarization  from  the  relative 
phases  of  its  components  recorded  along  three  orthogonal  seismometer  axes,  Polarization  places  constraints 
on  the  sometimes  messy  job  of  identifying  phases  present  in  the  coda,  most  challenging  in  regional  data 
with  closely-spaced  arrivals  and  scattering  (Gal’perio,  1984).  For  example,  only  S  waves  can  be  circularly 
or  ellipiically  polarized,  and  Rayleigh  surface  waves  will  be  elliptlcaily  polarized  until  the  P  component 
attenuates,  when  it  becomes  vertically  linearly  polarized  (Gubbins,  72;  1990).  Also,  the  polarization  angle 
of  S  waves  relates  to  parameters  entering  the  description  of  S  radiation  patterns  (Kasahara,  45;  1981),  which 
could  tell  us  something  about  the  type  of  source.  Furthermore,  when  presented  as  a  flinction  of  frequency 
and  time  in  tandem  with  a  sonogram,  the  polarization  helps  to  visually  locate  patterns  in  energetic  arrivals 
due  to  such  factors  as  dispersion,  scattering,  and  noise,  and  to  spell  out  an  arrival’s  region  and  distribution 
of  influence  in  both  time  and  frequency. 

To  estimate  the  principal  polarization  vector  (for  each  time  and  frequency  point),  we  first  fonn  a 
matrix  S  from  the  three  seismometer  component  column  vectors  whose  elements  are  the  less  leaky  wavelet 
transforms  (t.e.,  from  the  column  vectors  Q  that  correspond  to  each  seismometer  component): 


S  = 


/  *0,1  *0,»  *0,8  \ 


'*1  -1,1  *L- l.S  *1-1,3' 


=  (Gi  G,  Gj) 


where  we  have  L  transforms  with  good  leakage  properties.  If  we  apply  the  singular  value  decomposition 
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(SVD)  to  the  matrix  S,  wo  can  approximate  the  principal  polarization  vector  for  certain  regions  of  the 
time-frequency  plane  (Lilly  and  Paric.  1995).  The  SVD  breaks  S  down  into  a  product  of  three  matrices,  the 
first  and  last  of  which  will  now  be  square  and  unitary: 

S  =  UEV»  (3.1) 

The  first  three  column  vectors,  u;,  uj,  and  u$  of  U  are  orthonormal  eigenvectors  of  SS',  the  first  three 
column  vectors  of  V  are  the  eigenvectors  of  5*  S,  and  the  first  three  singular  values  of  the  diagonal  matrix  E 
are  the  corresponding  eigenvalues  (Strang,  appendix;  1988).  U  is  L*L,  E  is  1*3,  and  V  is  3*3.  Because 
the  first  three  columns  of  U  are  eigenvectors  of  S  Sf ,  they  are,  in  effect,  spectral  estimates  along  the  wave’s 
principal  axes  of  motion.  Rotating  S  by  a  certain  matrix  R  gives  these  spectral  estimates  weighted  by  the 
eigenvalues: 

SR  =  (*iUi  <rjuj)  (3.2) 

But  this  matrix  R  is  none  other  than  V  itself  (since  V  is  unitary),  which  means  that  the  columns  of  V 
contain  the  direction  cosines  between  the  principal  axes  of  motion  (i.e.,  of  wave  disturbance,  not  necessarily 
of  propagation)  and  the  seismometer  axes.  Furthermore,  if  a\  »  <rj,  0-3,  then  most  of  the  motion  occurs 
along  ffiui,  in  which  case  vi  denotes  the  principal  polarization  vector  (Lilly  and  Park,  1995).  These 
polarization  estimates  are  good,  however,  only  for  certain  regions  of  the  time-frequency  plane,  since  they 
are  predicated  upon  a  relatively  large  first  singular  value,  a  condition  satisfied  only  in  certain  regions.  It  is 
possible  to  assess  tbeir  validity  by  calculating  and  plotting  confidence  estimates  of  the  first  singular  value 
based  on  Monte  Carlo  simulations  of  the  SVD  (Lilly  and  Park,  1995). 

3.2  Results  for  synthetic  data 

The  accuracy  and  power  of  the  polarization  calculations  is  tested  on  three  synthetic  cases  whose  ex¬ 
pected  polarizations  are  known:  a  two-component,  completely  in-phase  sinusoidal  pulse,  a  two-component 
ninety  degree  out-of-phase  sinusoid,  and  a  single-component  chirp  signal,  with  a  small  dose  (3  %  of  the 
maximumsignal)  of  random  noise  added  in  all  cases  (see  Figures  3.1  through  3.8).  The  noise  is  added  to 
strengthen  the  test,  for  we  expect  to  see  random  polarizations  over  lime  and  frequency  bandwidths  containing 
noise  and  coherent  polarization  over  time  and  frequency  regions  containing  predominantly  coherent  signal. 
Clearly  in  the  first  case  (Figures  3.1  through  3.3),  we  expect,  over  the  lime  and  frequency  range  of  the 
energy,  a  linearly-polarized  signal  oriented  at  a  45  degree  angle  to  both  orthogonal  components,  and  random 
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In-phase  two-component  synthetic  sinusoidal  pulses,  noise  added 
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Figure  3.1  The  two-component  signal  corresponding  to  a  noisy  sinusoidal  poise  linearly* 
polarized  in  the  transverse-radial  plane  along  with  pure  noise  along  the  third,  vertical 
direction. 


polarization  due  to  noise  elsewhere.  In  the  second  case  (Figures  3.4  through  3.6)  we  expect  a  circularly- 
polarized  signal  over  the  time  and  frequency  intervals  coincident  with  coherent  energy,  and  again  random 
polarization  elsewhere.  For  the  chirp  input  (Figures  3.7  through  3.8)  we  expect  a  dispersive  spectrum  with 
linear  polarization  over  regions  of  significant  energy,  and  randomly-polarized  noise  everywhere  else.  This 
is  indeed  what  we  observe  in  all  cases  as  Figures  3.7  and  3.8  il  lustrate.  The  grid  of  lines  or  circles  appearing 
over  small  time  and  frequency  bands  and  superimposed  upon  the  spectrum  represent  the  signal  polarization 
for  the  given  time  and  frequency  band.  The  component  of  polarization  that  lies  parallel  to  the  horizontal 
axis  represents  the  radial  component,  and  that  parallel  to  the  vertical  axis  represents  either  the  transverse  or 
the  vertical  signal,  as  noted  in  each  case.  The  polarization  locks  right  onto  the  expected  orientations  over 
energetic  times  and  frequencies  and  remains  random  over  regions  of  pure  noise. 
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Linearly  polarized  (l.e„  in  phase)  two-component  synthetic  sinusoidal  pulse  with  noise  added 
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time,  s;  polarization  is  vertical-radial 

Figure  3.2  The  vertical-radial  polarization  superimposed  on  the  radial  spectrum  of  the 
pulse  in  the  last  figure;  note  that  the  polarization  is  denoted  by  the  superposed  magenta 
curves,  with  the  radial  component  always  plotted  along  the  horizontal  axis;  the  spectrum 
amplitude  corresponds  to  the  color  in  the  time-frequency  plane. 


3.3  Identification  of  surface  wave  phases  using  polarization 

Figures  3.9  and  3.10  show  the  .OS  to  1  hi  sonogram  spectrum  and  polarization  (indicated  by  the 
overlying  purple  circles,  lines,  and  dots)  of  a  ripple-fired  quany  blast  at  the  station  LBOH  in  Wyoming. 
Figure  3.9  displays  the  vertical-radial  polarization  (radial  component  along  the  horizontal  axis,  vertical  along 
the  vertical  axis)  overlying  the  vertical  spectrum,  while  Figure  3.  lOdisplays  the  transverse-radial  polarization 
(now  the  transverse  component  along  the  vertical  axis)  with  the  transverse  spectrum.  These  results  provide  a 
good  example  of  the  power  of  the  polarization  calculation  to  determine  possible  identifications  of  energetic 
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Linearly  polarized  (l.e.,  In-phase)  two-component  synthetic  sinusoidal  pulse  with  noise  added 
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Figure  3.S  The  transverse-radial  polarization  of  the  linearly-polarized  pulse  over  the 
radial  (or,  equivalently,  the  transverse)  spectrum;  random  noise  corresponds  to  random 
polarization  where  the  siganl  is  insignificant. 

components.  The  ambiguity  present  in  the  arrival  times  of  different  phases,  complicated  by  dispersion, 
would  make  it  difficult  to  identify  them.  It  also  provides  a  much  more  accessible  means  of  identifying 
phases  and  their  group  velocities  than  petforming  the  commonly-used  Multiple  Filter  Analysis  proposed 
by  Dzlewonski,  et  al.  (1969).  Without  the  polarization  information,  one  might  be  tempted  to  ascribe  the 
earliest-arriving  energy  to  SH  surface  waves,  since  such  surface  energy  typically  travels  faster  than  that  of 
P-S  V  surface  waves.  But  in  fact  we  see  with  the  aid  of  the  polarization  curves  that  the  SH  energy  arrives 
simultaneously  with  that  of  P-SV,  and,  though  generally  sandwiched  at  frequencies  below  it,  overlaps  it  to 
some  extent,  lb  distinguish  the  transverse  and  vertical  components  of  spectra  overlapping  one  another  In  the 
time-frequency  plane  would  be  difficult  from  the  spectra  alone.  Moreover,  though  multi-component  spectra 
may  pinpoint  regions  of,  say,  simultaneous  vertical  and  radial  energy,  they  cannot  tell  us  whether  the  energy 
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Synthetic  pi/2  out-of-phase  two-component  sinusoidal  time  series,  noise  added 
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Figure  3.4  The  two-component  lime  series  of  a  synthetic  noisy  circularly  polarized 
sinusoidal  signal  along  the  radial  and  transverse  directions  with  pure  noise  along  the  third, 
vertical  direction. 

is  ellipticaily  or  linearly  polarized,  and,  in  either  case,  at  which  angle.  These  are  also  very  important  pieces 
of  information  for  phase  identification.  We  would,  for  example,  expect  a  Rayleigh  wave  to  be  ellipticaily 
polarized,  and  we  might  find,  as  another  example,  that  a  particular  polarization  angle  conveys  information 
about  the  direction  of  force  at  the  source,  or,  alternatively  if  this  direction  is  known,  the  source-to-recelver 
azimuth  (see  Gal’perin,  1984,  and  Kasabara,  1981). 
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Synthetic  circularly-polarized  (i.e.,  pl/2  out-ot-phase)  two-component  signal  with  noise  added 


time,  s;  transverse-radial  polarization 


Figure  3.5  The  transverse-radial  polarization  and  the  radial  spectrum  of  the  circularly 
polarized  signal. 
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Figure  3.6  The  vertical-radial  polarization  and  radial  spectrum  of  the  circularly  polarized 
signal. 
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One-component  chirp  synthetic  with  noise  on  all  components 
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Figure  3.8  The  transverse  (or  vcriica])-radial  polarization  and  radial  spectrum  of  ihc  noisy 
synthetic  chirp  signal. 
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July  19  1996,  Wyoming,  Iboh;  vertical-radial  polarization,  vertical  spectrum 
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Figure  3.9  'Hie  venical-radial  polarization  and  vertical  spectrum  of  the  ripple-fired  event 
recorded  in  Wyoming  by  LB  OH. 
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July  19  1996,  Wyoming,  Iboh;  transverse-radial  polarization,  transverse  spectrum 


